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Dispersion in time and space: 

what propagating optical pulses in time (& not space) forces us to confront 

P. Kinslefl 

Blackett Laboratory, Imperial College London, Prince Consort Road, London SW7 2AZ, United Kingdom. 

(Dated: Friday 20 th February, 2015) 

I derive a temporally propagated uni-directional optical pulse equation valid in the few cycle limit. Temporal 
propagation is advantageous because it naturally preserves causality, unlike the competing spatially propagated 
models. The approach generates exact coupled bi-directional equations, which can be efficiently approximated 
down to a uni-directional form in cases where an optical pulse changes little over one optical cycle. It also also 
allows a direct term-to-term comparison of an exact bi-directional theory with an approximate uni-directional 
theory. Notably, temporal propagation handles dispersion in a different way, and this difference serves to 
highlight existing approximations inherent in spatially propagated treatments of dispersion. Accordingly, 
I emphasise the need for future work in clarifying the limitations of the dispersion conversion required by 
these types of approaches; since the only alternative in the few cycle limit may be to resort to the much more 
computationally intensive full Maxwell equation solvers. 

NOTE: despite the implication of the claimed “substantial text overlap” flagged up by arXiv, the phys¬ 
ical meaning of this work differs from that of arXiv:0810.5689 - it merely uses the same technique for 
decomposing waves into directional pieces, and as a consequence can have similar phrasing even when the 
physics is distinct. 


I. INTRODUCTION 

The importance of having simple and robust methods for 
the propagation of optical pulses has undergone increasing at¬ 
tention in recent years. This is due to the multitude of appli¬ 
cations 111] in which increasing short pulses are relied upon to 
act as a kind of strobe-lamp to image ultrafast processes HH, 
or even designer electric field profiles M to excite specific 
atomic or molecular responses. Alternatively, strong nonlin¬ 
earity is used to construct equally wide-band but also tempo¬ 
rally extended pulses - i.e. (white light) supercontinnua |7 l >| 
- or even come full circle and use the strong nonlinearity to 
generate sub-structure that is again temporally confined, as in 
optical rogue waves or even the temporally and spatially 
localized filamentation processes 11 1 [Hi¬ 
lt is clear, therefore, that progress towards diminishing 
pulse durations as well as their increasing spectral band- 
widths, and higher pulse intensities are all either stretching 
existing pulse propagation models to their limits, or breaking 
them. In such regimes, we need to be sure that our numeri¬ 
cal models still work, and have a clear idea of what has been 
neglected, and what the side-effects of those approximations 
are. In this paper I use a modern directional approach to show 
how a straightforward and relatively simple derivation allows 
a side-by-side comparison of exact and approximate propaga¬ 
tion equations, whilst still providing the numerical and ana¬ 
lytical convenience of a first-order wave equation. 

However, unlike the typical approach in nonlinar optics, I 
consider propagation in time rather than along some chosen 
spatial axis CIEI. The resulting evolution equations for 
the wave profiles in space no longer give us access to the 
frequency spectrum 11611 . but instead we have the wavevec- 
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tor spectrum, as well as a causally appropriate time evolution 
jl7|. This requires an alternate approach to the temporal re¬ 
sponse of the propagation medium, which highlights the ex¬ 
isting and somewhat poorly characterised approximations in¬ 
herent in spatially propagated treatments of dispersion. 

I start in section [TT] by deriving a customized second or¬ 
der wave equation from Maxwell’s equations, and reorga¬ 
nize it to define the material properties appropriately and set 
up the factorization stage. In section [TTT] I use the method 
of factorization that allows us to construct an explicitly bi¬ 
directional model, and which is then reduced to the popular 
uni-directional limit in section||V] Section lTY Dl discusses typ¬ 
ical modifications that can be applied to the equations given 
in sections uni and [TV] in order to and simplify them appro¬ 
priately and compare them to existing models. Having de¬ 
rived this main result - the propagation equations and their 
various specialisations - I turn in [3 to the consequences of 
temporal propagation and the new perspective it provides on 
the handling of dispersive effects in this and standard (spa¬ 
tially propagated) models. This is then followed in sectionlYTI 
by a specific comparison between the ordinary spatially prop¬ 
agated nonlinear Schrodinger (NLS) equation - perhaps the 
most widely investigated equation in nonlinear optics - and 
its temporally propagated counterpart. Finally, in section [VTT1 
I present my conclusions. 


II. SECOND ORDER WAVE EQUATION 

Most optical pulse problems consider a uniform and source 
free dielectric medium. In such cases a good starting point is 
the second order wave equation for the electric field, which re¬ 
sults from the substitution of the V x H = r), D |- J Maxwell’s 
equation into V x E = — d,B (see e.g. [Ts']). Here, however, 
I want to follow the displacement field D (not E) since it nat¬ 
urally occurs in conjuction with a time derivative - likewise. 
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the other important field is B rather than H. Thus we rexpress 
the curl Maxwell’s equations as 

Hod, D = ji 0 V xH-^ioJ = V x B — n 0 V x M - jU 0 J- 

( 2 . 1 ) 

£o(9fB = — £qV x E — EqK = —' V x D + V x P- e 0 K. 

( 2 . 2 ) 

Note that for both of these equations, we will also (as usual) 
want to specify the material response, whether dielectric or 
magnetic, in terms of a dielectric polarization P = D CoE 
and a magnetization M = B /ioH. Either of these equations 
given above might be either in the spatial domain (with ar¬ 
gument r = {r x . r v , r-)), or the wavevector domain (with argu¬ 
ment k = (k x ,k y ,k z )). 

Now we define the total magnetization M in terms of a 
component Mb linearly dependent on B, and another com¬ 
ponent Mb containing the rest. As a result, the excess polar¬ 
ization Mj, = Mb, which will typically include any nonlinear¬ 
ity, is also a function of B and not H. In traditional H-based 
spatially-propagated approaches, the temporal response (“dis¬ 
persion”) can be incorporated naturally, here it cannot. We 
have that 

M(r)=M L (r)+MB(r) = a M (r)B(r) +M B (r), (2.3) 

which in the spatial frequency (wavevector k) domain is 

M(k)=M L (k)+MB(k) =a M (k)*B(k)+MB(k), (2.4) 

since after a Fourier transform from r into k, the products be¬ 
come spatial convolutions (“★”); which emphasises that any¬ 
thing wavevector dependent is necessarily a non-local con¬ 
cept. Note that this is the converse of the usual situation where 
a temporal response is described using time-domain convolu¬ 
tions, and these become products in the frequency domain. 

Now we define the total polarization P in terms of a com¬ 
ponent Pb linearly dependent on D, and another component 
Pd containing the rest. As a result, the excess polarization 
P e = Pd, which will typically include any nonlinearity, is also 
a function of D. We have that 

P(r)=P L (r)+P D (r) = Cfe(r)D(r) +P D (r), (2.5) 

which in the spatial frequency (wavevector k) domain is 

P(k) = a E (k)*D(k) +P D (k). (2.6) 

With these definitions for <9,D, <9 ; B. and for M and P set 
up, we can now proceed with the derivation of a second order 
wave equation for D - 

IM)d,V = V x B - £i 0 V x M L - jU 0 V x Mb - MoJ (2.7) 

= VxB-Vxa ( ,B-/ioVxMB-|ioJ (2.8) 

= V x B - V xB - (Vau) x B x M s - jUqJ 

(2.9) 

= Vx [l-a M ]B-ji 0 VxMB-/ioJ- (2.10) 


Now take the time derivative of this equation, use the fact that 
afj is independent of time, define dp = 1 — dp and a E = 1 — 
a £ ; and then substitute for <9,B, so that 

£oAto<5, 2 D = £ 0 V x a^9 ,B - £ 0 ftoV x <9 ,Mb - EoMo^/J (2.11) 
c~ 2 d 2 D = -V xcif, [VxD-Vx P + £oK] 

- c“ 2 V x (9 ,Mb — c~ 2 d,J (2.12) 

d 2 D = —c 2 V xa^ [VxD-Vx a e D - VxP D + £oK] 

— V x <9 ,Mb — d,J (2.13) 

= -c 2 V xdfi [Vx(l-a £ )D-VxP D + £oK] 

- V x <9 ,Mb - d,J (2.14) 

= —c 2 V x dpV x a f D + c 2 V x x Pd 

- c 2 V x dfiE oK - V x <9 ,Mb — d ,J 

(2.15) 

Remember that in this picture, all field properties are spec¬ 
ified as functions of our chosen primary field D; notably, the 
magnetic field B which is usually the main argument for Mb 
must be determined from the calculated D using the relevant 
Maxwell’s equation (i.e. eqn. ( 12.21) ). 


A. Simple case: isotropic & homogeneous 

Since trying to treat all the details of eqn. (12.151 ) correctly 
leads to excessively complicated expressions, I will first sim¬ 
plify it down into the case where (a) the material’s reference 
properties dp, a £ do not vary in space, and (b) the effective 
monopole current is zero (K = 0). Since (magnetic) mono- 
pies do not exist, and all non reference properties can still 
be encoded in Pd, Mb, or J, this simplification does not (of 
itself) impose any approximation. However, the first simpli¬ 
fication means that our reference properties cannot incorpo¬ 
rate any exact knowledge about the spatial structure of the 
propagation medium 1 . This restriction on obtaining the best 
possible match between the reference behaviour and the ex¬ 
act behaviour of the propagation medium will typically affect 
later approximations, where we assume deviations from the 
reference behaviour are small. One additional advantage of 
selecting constant reference parameters a,, and a £ is that we 
can replace them with matrices to allow for e.g. birefringence 
and cross-polarization couplings. 

The simplified second order wave equation for changes in 
the displacement field D(r) as it propagates forward in time 


1 This differs from the situation present in the complementary spatially prop¬ 

agated approach |l4], where the medium’s linear temporal response typi¬ 
cally can be subsumed into the reference behaviour. Nevertheless, in nei¬ 
ther case can the spatial properties be incorporated exactly into the ref¬ 
erence behaviour; although approximations - notably that of a so-called 
"spatial dispersion” - can be made. 


2 







^0 


D20WE 


Dr.Paul.Kinsler@physics.org 
Time propagated optics and dispersion http://www.kinsler.org/physics/ 


is 

<9 2 D = —c 2 a ll a £ 'V x V x D +c 2 a^V x V x Pd 

— V X d t Ms — <3rJ (2.1) 

= +c 2 a M fl £ V 2 D — • D 

— c 2 fl M V 2 P 0 +c 2 a fl VV -P 0 

- V x d, M b - d,i. (2.2) 

Note that the retention of the V • D and V • Pd terms allows us 
to include charge effects, as is needed when discussing high- 
power nonlinear situations such as filamentation 0. 

In the spatial Fourier “wavevector” regime, we replace V —>■ 
;k, so that k 2 = k ■ k, and reorganize slightly. The equation 
for changes in the displacement field D(k) as it propagates 
forward in time is 

<3, 2 D — c 2 £fya E ZrD = —c 2 ajjk 2 Po — d t J — tk x 

— c 2 an a £ kk • D -f c 2 a /1 kk ■ Pd (2.3) 
<3, 2 D — D 2 D = +Ora £ 'Pd — <9,J — ;k x <9 ,Mb 

- £l 2 k 2 kk D + Qrk~ 2 a ~ 1 kk P D , 

(2.4) 

where we have defined a reference frequency Q for the prop¬ 
agation, based on the reference material properties, i.e. 

f!. 2 (k) = c 2 a M a £ k 2 . (2.5) 

We could - if we believed we understood the spatial proper¬ 
ties of the propagation medium sufficiently well - simply alter 
this definition by incorporating an alternate k dependence that 
mimics an assumed spatial dispersion. However, it needs to 
be understood that such a step is distinctly ad hoc, since it 
lacks the convolutions over k necessary for any accurate rep¬ 
resentation of the material structure. But, that said, such an 
approach may still have considerable practical utility; and in¬ 
deed is functionally equivalent to how other methods handle 
spatial dispersion, as a non-local (and hence, strictly speak¬ 
ing, a non-causal) process: e.g. those that import a waveg¬ 
uide dispersion calculated from the (spatial) modes given by 
the guide’s transverse cross-section. To allow for this possi¬ 
bility, in the following I will replace the physically rigorous 
magnitude k argument for Q, with k, its vector counterpart 2 . 
This allows for not only more complicated assumed spatial 
dispersions, but also additional orientation dependence. How¬ 
ever, these benefits (only) arise as a result of approximation. 
Strictly speaking, both on the grounds covered above, and 
ones taking stricter view of causality 0, any notion of spa¬ 
tial dispersion is incompatible with locality. 

Unlike for the spatially propagated scheme in (l3l, we have 
not included temporal response of the propagation medium in 


2 In some cases, depending on the approximate model chosen, this new t2(k) 
definition may also need to be convolved with other parts of the expression, 
e.g. Qr(k) Po(k) —> £2 2 (k) *Po(k), although such an additional complica¬ 
tion is likely to prevent us obtaining much conceptual or mathematical ben¬ 
efit. Consequently, in what follows I do not allow for such convolutions. 


our definition of the directional fields, since it is incompatible 
with our temporal propagation. To include this we must model 
it explicitly as is done in Maxwell solvers such as FDTD |2(t 
121]; or approximate it as if it were a spatial effect (see e.g. 
[2211, or the discussions in the later section IV1 


III. FACTORIZATION 


I now factorize the wave equation lfT4Ll23h . a process which 
neatly avoids some of the approximations necessary in tradi¬ 
tional approaches. Factorization takes its name from the fact 
that the LHS of eqn. (12.2b or (12.4b is a simple difference of 
squares which might be factorized, indeed this is what was 
done in a somewhat ad hoc fashion by Blow and Wood in 
1989 m. Since the factors are just of the form (<9, =F ik), each 
by itself looks like a forward (or backward) directed wave 
equation. It therefore allows us to define a pair of counter- 
propagating Greens functions, and so divide the second order 
wave equation into a bi-directional pair of coupled first order 
wave equations. In addition, it allows us to directly compare 
the exact bi-directional and approximate uni-directional theo¬ 
ries term for term, whereas in other approaches the backward 
parts simply vanish and are not directly available for compar¬ 
ison. A short summary of the factorization method is given in 
appendixlAl 

To re-iterate an important point - the choice of a £l a^ and 
therefore £2(k) in eqn. (12.4b defines the specific Greens func¬ 
tions used; it therefore also defines the basis upon which we 
will then propagate the displacement field D. 


A. Bi-directional wave equations 


A pair of bi-directional wave equations suggests similarly 
bi-directional fields, so I split the electric displacement field 
into forward (D+) and backward an directed parts, with 
D = D + +D~. In the following equations, remember that Pd, 
Mb, etc all depend on the full field D, and not only one par¬ 
tial field (e.g. only D + ). This is an important point since we 
see they drive both forward and backward evolution equations 
equally. 

The coupled bi-directional first order wave equations for 
the directed fields D + are propagated forward in time whilst 
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being evolved forward and backward in space are 


d,D±(k) = ±^(k)D ± (k) ± ' Po(DT;k) 

T 2tt(k) kx “'’ (D -' ;k| 

T/“l (k .lkk D(k) ± in ' ( !‘ ) °t‘ kk-P J (D,»;k) 


2£2(k)k 2 


<9,D ± (k) =±;f2(k)D ± (k) ± 


;£2(k) 

2 a £ 


2£2(k )k 2 


Pz>(D,/;k) 


(3.1) 


^mr) ^m) kxMa(D - nk) 

T «“. D(k)± ^i|.p D(D , :k) 


2 k 2 


(3.2) 


Here I have replaced the partial time derivatives on the RHS 
with over-dots, this being to emphasise that if our expression 
is to remain strictly causal 10. then our models for current 
J and magnetization M® must return their time-derivatives as 
explicit functions of known quantities (here, typically D or t). 

Note that factorized equations can be rebuilt into a single 
second-order equation by taking the sum and difference, then 
substituting one into another with the assistance of a derivative 
w.r.t. t (see section IV.B of ll23ln . 

A remaining clarification required is to understand the 
meaning of “forward” and “backward” in the context of three 
dimensional space - this being less straightforward than in the 
spatially propagated case, where one spatial axis is selected 3 , 
and forward and backward refers to time t. Here, we should 
understand “forward” and “backward” as being with respect 
to a particular choice of k, and those individual vector com¬ 
ponents, that we are interested in. 


B. Propagation, evolution, and directed fields 

Note that since our solutions of the wave equations enforce 
propagation towards later times f, the fields D ± (t) are directed 
forwards and backwards in space; these fields then evolve for¬ 
wards and/or backwards in space as 1 increases. Note that I 
use this terminology (propagated, directed, evolved) through¬ 
out this paper to mean these three specific things. 

When examining the wave equation eqn. ( 13.21 which 
evolves the directed fields D ± as they propagate forward in 
t, we see that the RHS has two types of terms: which I label 
the “reference” and “residual” parts |25{]. 

Reference evolution (or “underlying evolution") is that 
given by ±z£2(k)D ± term, and is determined by our chosen 
a £ and cq,, and any additional ad hoc refinements. By itself, 
it would describe an ordinary oscilliatory evolution where the 


3 Typically a cartesian axis, although other choices are possible 0. 


field oscillations in I) (r) would move forward (+) or back¬ 
wards (—) in space. This is analogous to the choice of ref¬ 
erence when constructing directional fields [26], but here is 
done for temporal rather than spatial propagation [27]. 

Residual evolution accounts for the discrepancy between 
the true evolution and the underlying/reference evolution, and 
is every part of the material response not included in a e ,a^ 
(and hence £2(k)). i.e. it is all the remaining terms on the 
RHS of eqn. (13.2I >. These typically include any non-linear 
polarization, or spatially or orientationally dependent linear 
terms; they are analogous to the correction terms used in di¬ 
rectional fields models. In the language used by Ferrando et 
al. [23], these residuals are “source” terms. Although we 
might hope they will be only a weak perturbation, so that we 
might make the (desirable) uni-directional approximation dis¬ 
cussed later, the factorization procedure is nevertheless valid 
for any strength. 


C. Reference evolution: choice of £2 and the resulting D + 

I now examine how the choice of £2 affects the relative sizes 
of the forward and backward directed D + and I) . To do this I 
consider the simple example of a medium for which the field is 
known to propagate with frequency ft); but for demonstration 
purposes we choose a reference evolution determined by a fre¬ 
quency £2 that is different from ft). E.g., for a linear isotropic 
medium we could exactly define ft) 2 = £2 2 + A 2 ; but in gen¬ 
eral we would just have some residual (source) term Q. This 
means that our definitions of forward and backward directed 
fields do not exactly correspond to what the wave equation 
will actually evolve forward and backward as we propagate 
towards larger t. 

The second order wave equation is (d 2 + £2 2 )D = —Q, 
which in the linear case has Q = A 2 D, so that (<9, 2 + ar)D = 0. 
The factorization in terms of £2 is then 

<9 r D ± = ±;£2D ± ±—. (3.1) 

2 £2 

Now if we select the case where our field D only evolves 
forward, we know that D = Do exp [zft)f]. Consequently />* 
must have matching oscillations: i.e. D ± = Djf exp[jft>t], even 
though D is directed backwards. Substituting these into eqn. 

ED gives 


Do 


£2- G) 
£2 + 0 ) 


D, 


(3.2) 


which specifies how much I) we need to combine with D + so 
that our pulse evolves forward; since the D will be dragged 
forward by its coupling to D + . This interdependent I) 1 be¬ 
haviour is generic - no matter what the origin of the discrep¬ 
ancy between £2 and the true evolution of the field (i.e. the 
residual/source terms such as mismatched linear behaviour, 
nonlinearity, etc): some non-zero backward directed field I) 
must exist but still evolve forwards with D + . Analogous 
behaviour can be seen in the directional fields approach of 
Kinsleretal. li26h . 
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Usually we hope that this residual D contribution is small 
enough so that it can be neglected. If we assume I) ~ 0, then 
we find that ft) ~ £2 + A 2 /2£2, which is just the expansion of 
CO = (£2 2 + A 2 ) 1 / 2 to first order in A 2 /£2 2 . Following this, we 
find that eqn. (13.2b then says that Dq ~ (A 2 /4£2 2 )Dq\ which 
has come full circle and provided us with the scale on which 
D can be considered negligible. Outside the restricted (lin¬ 
ear) case where we know A 2 , the true frequency <o might be 
difficult to determine, and in nonlinear propagation any local 
estimate of the frequency can change as the pulse propagates. 

There is a further important point to notice: if we choose 
£2 = £2(k) with a wave-vector dependence, then we see that 
the source-like terms inherit that dispersion; and this is almost 
inevitable, since typically £2 ~ ck. This means, for example, 
that even if we started with nonlinear polarization model with¬ 
out any k dependence, our factorized equations will have a k 
dependence on the nonlinear terms. This is the complemen¬ 
tary process to what happens in the spatially propagated direc¬ 
tional fields approach of Kinsler et al. |2(|, where choosing a 
dispersive reference has an equivalent effect on the correction 
terms. 


IV. UNI-DIRECTIONAL WAVE EQUATIONS 

Making only a single, well defined class of approximation 
I can now reduce the exact coupled bi-directional evolution 
of D down to a single uni-directional first order wave equa¬ 
tion. As with the complementary spatially propagated theory 
01,1 do not require a moving frame, a smooth envelope, or 
to assume inconvenient second order derivatives are somehow 
negligible, as was necessary in traditional treatments. The ap¬ 
proximation is that the residual terms are weak compared to 
the (reference) ±/£2D term - e.g. weak nonlinearity or ori¬ 
entation dependence. This enables me to assert that if I start 
with D~ = 0, then I) will remain negligible. In this context, 
“weak” means that no significant change in the backward field 
is generated in a distance shorter than one wavelength (“slow 
evolution”); and that small effects do not build up gradually 
over propagation times of many periods (“no accumulation”). 

Slow evolution is where the size of the residual terms is 
much smaller than that of the underlying linear evolution - 
i.e. smaller than £2D. This allows us to write down straightfor¬ 
ward inequalities which need to be satisfied. It is important to 
note the close relationship between these and a good choice of 
£2, as discussed in subsection lHICI If £2 is not a good enough 
match, there always be significant contributions from both for¬ 
ward and backward directed fields; and even if nothing ends 
up evolving backwards, an ignored backward directed field 
will result in (e.g.) miscalculated nonlinear effects, since the 
total field D = D + + D will be wrong. 

No accumulation occurs when the evolution of any back¬ 
ward directed field is dominated by its coupling via the 
residual terms to the forward directed field D + ; and not by its 
preferred underlying backward evolution. No accumulation 
means that forward evolving field components do not couple 
to field components that evolve backwards; this the typical 
behaviour since the phase mismatch between forward evolv¬ 


ing and backward evolving components is ~ 2£2; in essence 
it is comparable to the common rotating wave approximation 
(RWA). This rapid relative oscillation means that backward 
evolving components never accumulate, as each new addition 
will be out of phase with the previous one; it is not quite a “no 
reflection” approximation, but one that asserts that the many 
micro-reflections do not combine to produce something sig¬ 
nificant. An estimate of the terms required to break this ap¬ 
proximation are given in appcndixlB] generally speaking this 
is a much more robust approximation than the slow evolu¬ 
tion one. Of course, a periodic temporal modulation of the 
medium will give periodic residual terms, and these can be 
engineered to force phase matching. In spatially propagated 
analyses, where we talk of spatial (wavevector) phase match¬ 
ing and not this temporal (frequency) phase matching [28], 
this would be a periodicity based on a relatively small phase 
mismatch (see e.g. quasi phase matching in Boyd |2S}[); but 
might even go as far as matching the backward wave (see e.g. 

It is also important to note that the same small size of pertur¬ 
bation from the residual terms can accumulate on the forward 
evolving field components (or, indeed, the backward pertur¬ 
bation on the backward evolving field components). Although 
the magnitude of the residual terms acting on the forward and 
backward field evolution are identical, forward evolving com¬ 
ponents of the residuals can accumulate on the forward evolv¬ 
ing field because they are phase matched; whereas backward 
residuals are not, and rapidly average to zero. 

A. Residual terms 

The handling and criteria for different types of residual 
terms has been discussed already 1 14-1 for the spatially prop¬ 
agated case. While it would be straightforward to repeat that 
level of detail here, in fact the basic process is virtually iden¬ 
tical except for the different roles played by the t and r argu¬ 
ments. Accordingly I will consider only the most commonly 
considered residuals - the dielectric scalar and vector terms 
as well as currents. Those interested in magnetic or magneto¬ 
electric effects can adapt the process to their own systems in a 
similar manner. 

The relevant total polarization P d can then be decomposed 
into pieces before seeing how each might satisfy the slow evo¬ 
lution criteria. Thus we write 

Po(D,r,f) = a e (D,r,f)D(r,f) +U D (D,r,f) 

= a L (r,t) *D(r,/) + ajv(D,r,f) *D(r,f) 

+ U L (r,0*D(r,0+U;v(D,iy). (4.1) 

The part which is scalar in nature is represented by a e , and 
might contain (non-reference) linear parts and time response 
(ai), perhaps a convolution (“*”) over the past, or angle de¬ 
pendence; it might also contain nonlinear contributions on 
such as the third order Kerr nonlinearity with «\T) °c (D • D)D. 
The vector part U/j, if present, could be from e.g. a sec¬ 
ond order nonlinearity, which couples the ordinary and extra¬ 
ordinary field polarizations. Note that this description of the 
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material parameters does not restrict them in any way, and can 
e.g. include any order of nonlinearity. 

B. Slow evolution? 

Now I will treat each possible residual term in order, where 
the oppositely directed field is negligible: i.e., for I)', we have 
that I) ~ 0. where the scalar contains the linear response 
of the material that is both isotropic and lossless (or gain-less); 
since here it is a time-response function, it is convolved with 
the electric displacement field D. Note that the field vector D, 
and indeed the material parameter El are all functions of time t 
and space r = (x,y,z)', the polarization P e and its components 
0 e , V e are a functions of time t, space r, and the field D. 

Below I will refer to field components Z), and Df, where 
D + + D = D = (D x .Dy.D z ) and i £ {x,y,z}; also to wave 
vector components kj from k = (k x ,k y ,k z ), with kf = k 2 + k 2 . 
However, note that in the constraints below, that kj is also used 
as a substitute symbol to represent any one of k x , k y , or k-. 

Any corrections due to imperfect choice of Q will appear 
in the Pd polarization correction term. However, unlike in the 
spatially propagated case, there is no diffraction-like term - 
here, that is quite naturally part of the reference behaviour. 

Firstly, scalar polarization terms (j), which can be either linear 
(0l(&)) or nonlinear (0^(D ,k)). These might encode e.g. the 
time-response of the medium (and hence its dispersion), or 
be some nonlinearity response such as the third-order Kerr 
nonlinearity. Since these are locally correlated in space, in the 
wavevector picture they include convolutions over k, and so 
give us the criterion 


lQ.(j)-k 

1 D] 

f + D/ 

\/2 <j)* 

1 Df 

I/ 2 

iQ. \Df 


Df | 



In the linear case, <j> = <f>L is independent of D, so only the ma¬ 
terial parameters are constrained and the pulse properties play 
no role. In the nonlinear case, the criteria makes demands on 
the peak intensity of the pulse - but does not apply smooth¬ 
ness assumptions or bandwidth restrictions. 


with the same strength, so that it is far from obvious that we 
can set D v to zero, but still keep the D y without being incon¬ 
sistent. However, as described above, it is the phase matching 
which ensures that forward residuals accumulate, whilst the 
non-matched backward residuals are subject to the RWA, and 
become negligible: hence we can still rely on eqn. (14.21) . albeit 
under caution. 

Thirdly, we have the residual terms indicated by the presence 
of time-varying currents. These are 

|/i| <2Q. 2 \Df\. (4.3) 

Fourthly, we have the divergence terms kk ■ D, kk ■ Pd, which 
amount to charge density contributions either from D or Pd, 
with ap(k) = aj 'k ■ Pd and Or>{k) = k ■ D. Thus if the 
wavevector is oriented along a unit vector u, then 

utjFl |oj /2k <C iQ\Df\ 

aui <C 2k |.d4|. (4.4) 

is the relevant criterion. 

To summarize, these various criteria assert that modulations 
away from the reference evolution must be weak. Notably, 
weak nonlinearity is invariably guaranteed by material dam¬ 
age thresholds 1 3111 . However, and in contrast to the usually 
favoured (but strictly speaking less causal) spatially propa¬ 
gated methods, the lack of any rigorous way to incorporate 
dispersion into the reference frequency Q,(k) makes it harder 
to satisfy the criteria for the ordinary linear response. 

C. Uni-directional equation for D + 

In the case where all of the slow-evolution criteria listed 
above hold, we can be sure that the backward directed field 
D is negligible, and if the no-accumulation condition also 
holds, then neither will there be any non-negligible backward 
evolving contributions to the field. Consequently, we can be 
sure that an initially negligible I) remains so, and eqn. (13.2b 
becomes 


Secondly, linear and nonlinear vector terms from U. These 
will have a criterion broadly the same as the scalar cases in 
eqn. (ED, but with U replacing ©I). Thus for i £ {x,y,z}, we 
can write down constraints for each component Uj, which are 

ia\Ui\/2<^ia\Df\ => |f/,-|<2|Z)f|. (4.2) 

In the linear case, U = U/., and since U/. and D have some 
linear relationship, this criterion only constrains the material 
parameters contained in U/., not the field profile or amplitude. 
The nonlinear case is similar to the scalar one, and the peak 
pulse intensity is restricted; e.g. for a % { ' 2) medium, U.v ~ 

x (2) |d|. 

However, one complication of the vector cases is that a field 
consisting of only one field polarization component (e.g. Df) 
may induce a driving in the orthogonal (and initially zero) 
components (e.g. Df). Hence both Df fields will be driven 


<3,D+(k) = +zft(k)D + (k) 
j(f;k) 


z£2(k) 

2 a £ 


Pd(D+, f;k) 


2H(k) 2£2(k) 

iQ.{ k) k 


kxM B (D+,f;k) 


2 k 2 


k-D+(k) +^^-^-k-PD(D + T;k). 


(4.1) 

where now the polarization Pd, diffraction, and divergence are 
solely dependent on the forward directed field D + . 


D. Modifications 

Just as for the comparable spatially propagated derivation 
| ll we might also apply some of the strategies used in other 
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approaches to get a more tractable and/or simpler evolution 
equation. Unlike in “traditional” derivations 02143511 . none of 
these are required, but they nevertheless may be useful. 

Briefly, a co-moving frame might be added, using z! = 
Z + tVf. This is a simple linear process that causes no ex¬ 
tra complications; the leading RHS icD + term is replaced by 
t(c=F v)D ± . Note that setting c = v will freeze the phase ve¬ 
locity of the pulse, not the group velocity. Also, a carrier- 
envelope separation could be implemented using D(z) = 
A(z)exp[i(fi>if — kjz)] +A*(z)exp[— i(cO\t -kiz)] defining the 
envelope A (z) with respect to wave vector k\ and carrier fre¬ 
quency co i. Multiple envelopes centred at different wave vec¬ 
tors and carrier frequencies (kj, (Oj ) could also be used to split 
the field up. This is typically done when the field princi¬ 
pally comprises a number of distinct narrowband components 
(e.g. II29L 136f 137[1). The wave equation can then be separated 
into one equation for each piece, coupled by the appropriate 
wavevector-matched polarization terms (cf |[28*]). 

V. DISPERSION 

Consider a dielectric slab waveguide of thickness d in the 
perpendicular x direction, propagation in the z direction, and 
with core and cladding permittivities £\ and £q. In the CW 
limit, the spatial structure of the waveguide gives rise to fre¬ 
quency dependent properties that are usually called “spatial 
dispersion”. Even for this simple slab design, the dispersion 
relation has a non trivial form; in this case it is given by the 
solution to a transcendental equation [ 3j$]. Typically, this is 
written in a way implying we want to calculate k x and ft = k- 
from ft), i.e. for the transverse electric (TE) field modes we 
have 

T (, k x d ) = kf 1 y/ ft) 2 Mo (£1 - £ 2 ) - £?, (5.1) 

where p 2 (co) = co 2 /1q£\ - k 2 xl (5.2) 

and T(k x d) is either t'&n{k x d) or — cot (k x d). 

However, we can now see that this is an implicitly spatially 
propagated representation, since in time propagated systems 
we do not (cannot) know ft), but we do know k. Therefore eqn. 
(15.2b should be rewritten accordingly, using k- for p (since 
wavevector is no longer a propagation reference), and f2 for 
co (since frequency now is a propagation reference); thus 


picture, we can still find the same waveguide modes and their 
dispersions, but this time given as functions of k z , not ft). So 
for some n-th mode of the waveguide, as expanded about some 
central evolution wavelength k zni which corrresponds to a ref¬ 
erence frequency £ 2 n = £ 2 (k zn ), we might have 

Q(k z ) =a,+L 7m (k z - k zn ) m , (5.5) 

m 

where yi = c + y 1 has two contributions - that due to the vac¬ 
uum (just c), and the estimated modifications due to spatial 
effects (/i). 

We may also wish to add the effects of whatever tempo¬ 
ral response the material might have to our propagation cal¬ 
culation; but also decide that an approximate model where 
we modify the expression above in eqn. (15.51) is sufficent. 
This is the converse of the traditional procedure, where the 
waveguide dispersion is added to whatever /3 (o) is appropri¬ 
ate for the medium’s time response. Consequently we need 
only make the same assumptions: (a) that the two sources of 
dispersion are weak and so can simply be added, (b) that the 
temporal response is assumed to have negligible effect in the 
waveguide’s propagating modes, and (c) that the temporal dis¬ 
persion is the same at all points in the waveguide. 

For simplicity we consider the usual situation where our 
dispersion has already been characterised and represented as 
a truncated power series in frequency. In line with the rec¬ 
ommendations of [25], I exclude imaginary components (e.g. 
loss or gain terms) in either the wavevector reference ft or the 
frequency reference £2. Consequently, the reference wavevec¬ 
tor considered here has the simple form 

P(co) = fio + ic, (co - coo) + k 2 (co - coo ) 2 , (5.6) 

where Kj = c _1 + kj has two contributions - that due to the 
vacuum (just 1 /c), and the modification due to the temporal 
response (kj). 

I have chosen this form for p (co) because it is possible to 
calculate its inverse, although in general - if there are more 
terms in the expansion, or /j has some specific and non-trivial 
analytic form - the problem is more complicated. For this we 
again replace ft) with fi to denote the new role as a reference 
frequency, and replace p with k z . With £.2 n = ck zn , this gives 


c 2 n\k x ) [T 2 (k x d) + 1]. (5.3) 

where k 2 = c 1 kl 2 (k x ) — k\. (5.4) 


k2-(kf) — kl n ^ 


1 T \/l +4(k : - k z „) k 2 / 


(5.7) 


Here we see that the natural (reference) propagation frequency 
is defined solely by the transverse wavevector k x , as is the 
matching evolution wavevector k z . While apparently simpler 
than the traditional form, since kl is directly given by k x , we 
usually want to neglect any direct knowledge of the transverse 
properties. This means that what we really want is not kl(k x ), 
but kl(k z ), the inverse of the traditionally preferred p (ft)). 

Despite these complications, which are in any case rather 
similar to those we see in the traditional spatially propagated 


which we could expand into a power series for small k zn — k z . 
Since Kj includes the vacuum, and dispersion k 2 is assumed 
to be small, we should make the top sign choice since it gives 
the physically relevant small frequency corrections. The lin¬ 
ear term in the series has the coeffient J\j = 1/Kj and the 
quadratic term Jn = Kj/Kj\ 

We now combine the two dispersive effects - the spa¬ 
tial/geometric/waveguide ones based on parameters 7 , and 
converted temporal ones based on Kj, to give a total reference 
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frequency based on k- of 


&'(,k z ) — + 7i {k z — k zn ) + Ym \{K — k m )] m , 

m> 1 


_Kj_ 

2k 2 



— £l n + 



1+4 (k z k zn ) Ki / fC| 

(k, — k, n ) + ^y 2 — —j ^ 


(5.8) 

(k; — k zn ) . 

(5.9) 


Note that to avoid incorporating the vacuum contribution 
twice when summing the dispersive effects, I have replaced 
7 i with fi ■ 

In the complementary case, useful for spatial propagation, 
we can follow the converse - and indeed the usual procedure 
- to convert the CW effect of the waveguide’s spatial structure 
kl(k z ) into a frequency domain form 1138113911 . This then can 
be combined with temporal response of )3 (ft)) to get 

j3'(ft)) = A>+ + (to ©o) + (k 2 -¥^J (co-coo) 2 . 

(5.10) 


Discussion 

We might have hoped to find a less approximate compari¬ 
son between the spatially propagated and the temporally prop¬ 
agated approaches to dispersion what has just been discussed. 
Further, although qualitative comparisons can also be made, 
ideally we would like a specific and quantitative comparison 
as done for e.g. testing the unidirectional approximation & 
That this is not as simple as it might sound can be demon¬ 
strated as follows: write two mathematically identical but 
physically and notationally distinct propagation equations, 
i.e. 


d z E{co) = +iK(co)E{co) (5.11) 

d t D{z) = +ikl{z)D(z). (5.12) 

For simplicity, this comparison is done in the 1+1D limit and 
for unidirectional equations; the residual ft) or z behaviour has 
been merged with the reference behaviour to give single prop¬ 
agation parameters K(co) or £l(z). 

Superficially, the comparison between eqn. (15.1 11 1 and 
(15. 12b looks promising. However, the vacuum part of K(co) 
is fi)/c, but the vaccum part of £2(z) is in contrast a spatial 
derivative, i.e. cd z . Even if it were possible to match up 
the functional forms of the non-vacuum temporal response of 
K(co) — co/c with a spatial structure giving £2(z) — cd z , the 
vacuum parts would remain unmatched. Alternatively, we 
could compare the z-propagated ft)-domain to the f-propagated 
/.-domain. However, although in this case the reference be¬ 
haviours match, but the /.-domain model now has a convolu¬ 
tion that does not appear in the co-domain one. 

As a result, it is unclear how well approximating temporal 
dispersion as spatial properties works in practise; but then it 
is equally unclear how well approximating spatial properties 


as temporal dispersion works - even if this latter process is so 
ubiquitous in optics as to be effectively invisible. Removing 
either of the mismatches discussed above - either in reference 
behaviour or convolution - requires a narrowband (CW) limit, 
a process which necessarily obscures the parameter regime 
where interesting tests and comparisons can (and should) be 
made - i.e. the short pulse, large bandwidth limit. Indeed, the 
stringent nature of the CW approximation would defeat the 
entire point of the comparison. 

Nevertheless, it may be possible to design numerical sim¬ 
ulations containing structures and temporal responses whose 
scales are carefully graduated so that the mundane bandwidth 
issues can be kept separate from the interesting propagation 
ones. Such a design would, as desired, allow the propagation 
approximations (only) to be compared and contrasted by nu¬ 
merical simulation. However, this is a non-trival task, and one 
which I leave for later work. 


VI. PULSE PROPAGATION 

Third order nonlinearities are common in many materials, 
e.g. in the silica used to make optical fibres 118.1381.1391. Here 
I compare 1 + 1D propagation such a waveguide, in order to 
elucidate the differences between the standard spatially prop¬ 
agated picture and the temporally propagated picture. 

The assumptions made are those of transverse fields, weak 
dispersive corrections, and weakly nonlinear response; these 
all allow us to decouple the forward and backward wave equa¬ 
tions. This decoupling then allows us, without any extra ap¬ 
proximation, to reduce our description to one of forward only 
pulse propagation. The specific example chosen here is for an 
instantaneous cubic nonlinearity, but it is easily generalized to 
non-instantaneous cases or even other scalar nonlinearities. 


A. Spatially propagated NLSE 


The spatially propagated NLS for Ej (ft)) based on a direc¬ 
tional factorization [141 can be written 


d z E+ = +ifiE+ +Y J Kn.((O-CO 0 ) m E 

m 

ik} 


2/r 


X {3) E 2 (t)E+(t) 


(6.1) 


where +[...] is the Fourier transform that converts the nonlin¬ 
ear polarization term into its frequency domain form. Here we 
have chosen a fixed, non-dispersive reference wavevector /3, 
the temporal response of the propagation medium is encoded 
in the coefficients Kj. For the temporal (time response) of the 
material, this is an in-principle exact re-representation of that 
response as a Taylor series expansion; although in practise the 
series is usually truncated after a few terms and the propagat¬ 
ing fields restricted to within a bandwidth where only those 
terms are significant. 

The important point to remember is that in this spatially 
propagated picture we do know the full time history, at least 
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in a computational sense - i.e. as we compute a solution to the 
propagation equation. This “computational past” [16], com¬ 
prising all previously computed values of the field properties 
and material state means that a full frequency spectrum and re¬ 
sponse is known. As a result, the time response of the material 
can be directly encoded in an accurate and useful dispersion 
relation, and then be approximated as a Taylor series in ft). 

Since this propagation is assumed to be taking place in a 
waveguide, we also have its geometric - wavelength sensitive 
- properties to consider. However, as discussed in the previ¬ 
ous section we typically just re-represent those as if they were 
instead generated by a temporal response (see e.g. 1138113911 ) 
Thus the coefficients Kj usually combine both the effect of 
time response and geometric properties into a single total dis¬ 
persion; this is the very basis of dispersion compensation in 
optical fibres, since the material (temporal) dispersion is off¬ 
set by the waveguide design (spatial dispersion) [38]. 


B. Temporally propagated NLSE 

Using the approach derived and presented in this paper, we 
can develop a time propagated counterpart to eqn. (16.11 ). In 
a purely mathematical sense, it is, apart from notation, iden¬ 
tical to eqn. (IQ ). However, in physical terms, the change 
is not trivial, since the physical meaning and boundary condi¬ 
tions differ significantly. Starting from eqn. (ED we can get 
a propagation equation in wavevector space for I) ■ (k) which 
reads 


d t D+ = +&D+ +£y m . (k z - k z0 ) m Dl 


iD. 

2 a? 


X^D 1 2 3 x (z)D+( z ) 


( 6 . 1 ) 


where 3 [•••] is the Fourier transform that converts the spatially 
localized nonlinear polarization into its wavevector domain 
form. Here, all dispersive behaviour has been combined into 
the coefficients Jj, which mimic the common notion of spatial 
dispersion. Such terms can arise directly from a Taylor ex¬ 
pansion in k of a relevant CW dispersion relation about some 
suitable reference wavevector; however, although we can use 
the expansion to represent the spatial structure of the propaga¬ 
tion medium (e.g. a waveguide), the conversion from spatial 
structure to k expansion is only relatively straightforward in 
the CW limit. 

If the propagation medium also has a temporal response 
(i.e. temporal dispersion), and if we do not want to model 


it explicitly, then we can use the approach described in the 
previous section, where the temporal response of the material 
was (approximately) represented as if it were instead due to 
spatial properties. 

Finally, note that in ordinary temporally propagated sim¬ 
ulations, the computational past and the physical (temporal) 
past are identical; this is not the case for spatially propagated 
simulations, as mentioned above. 

VII. CONCLUSION 


I have derived a general first order wave equation for uni¬ 
directional pulse propagation in time that allows for arbitrary 
dielectric polarization, diffraction, and free electric charge and 
currents. This propagation equation has a utility all of its 
own, allowing efficient unidirectional propagation in a strictly 
causal manner (i.e. time directed). However, it requires com¬ 
parable approximations to those used in deriving spatially 
propagated unidirectional propagation equations, although the 
specific tradeoffs are different. 

The propagation equation was derived by first factoriz¬ 
ing the second order wave equation for D into an exact bi¬ 
directional model, it applies the same type of well defined 
type of approximation to all non-trivial effects (e.g. nonlin¬ 
earity, diffraction), and so reduces the propagation equations 
to a first order uni-directional wave equation. One nice fea¬ 
ture of this factorization method is that it provides for a term- 
to-term comparison of the exact bi-directional theory with its 
approximate uni-directional counterpart. 

In addition to its use where the causal simulation of a prop¬ 
agation problem is of strict concern, it also illuminates the 
rather shadowy process of handling material response and spa¬ 
tial structure in waveguides by constructing a single “com¬ 
bined dispersion” parameter. As discussed in section [V] we 
can now see that this simple addition of temporal and spatial 
dispersions in an ad hoc manner stands on rather poor foun¬ 
dations, particularly in the few-cycle pulse limit. Further, it 
is hard to see any straightforward way of testing the limita¬ 
tions of this common procedure, which poses an interesting 
challenge for the future. That said, in most applications, we 
expect that a combined dispersion will remain a perfectly ad¬ 
equate approximation, since the pulses are rarely so short, and 
it seems reasonable to assume that the character of the pulse 
evolution should remain similar, even if details might differ. 
Lastly, the lessons learnt from this investigation could also 
be applied to other types of wave propagation, most notably 
acoustic systems (see e.g. [|271). 
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Appendix A: Factorizing 


First, consider a forward evolving field so D = Z>oexp(/ft)f), 
and therefore 


Here I present a simple overview of the mathematics of the 
factorization procedure, since full details can be found in [23;]. 
In the calculations below, I alter the physics but not the math¬ 
ematical process (of I 23 IP to instead transform into frequency 
space, where the /-derivative d t is converted to /ft). Also, we 
have that Or = c 2 k 2 /n 2 , and the unspecified residual term is 
denoted Q. The second order wave equation can then be writ¬ 
ten 


[d 2 + Q 2 ]D = -Q 
[-ft r + Qr]D=-Q 


-1 

2Q 


1 

ft) T £2 


(A.0.1) 
(A.0.2) 


1 

(ft) - £2) (ft) + £2) 

(A.0.3) 


1 

co-Q 


Q- 


(A.0.4) 


Now (ft) — £2) _1 is a forward-like (Green’s function) propaga¬ 
tor for the field, but note that in my terminology, it evolves 
the field. The complementary backward-like propagator is 
(ft) -flT)^ 1 . As already described in the main text, we now 
write D = D + + D~, and split the two sides up to get 


zr 


D~ = 


D ± = 


-1 

2 p 

±1 


1 


1 


ft) + £2 co — Q 
1 


Q 


[fi)Tfl]D =± 


2Q ft>q=£2 
1 1 


Q 


2Q ft>q=£2 


Q 


icoD ± = ±iQD + ± — Q. 

2Q^ 


(A.0.5) 
(A.0.6) 
(A.0.7) 
(A.0.8) 


Finally, we transform the frequency space /ft) terms back into 
normal time to give t derivatives, resulting in the final form 


=$ D o, (B-0.1) 

where as noted ft) can be difficult to determine, and may even 
change dynamically; here we can assume it corresponds to 
the propagation wave vector that would be seen at if all the 
conditions holding at a chosen position also held everywhere 
else. On this basis, we can (might) even (try to) define ft) = 
ft)(f), where by analogy to the linear case we might assert that 
ft) 2 (f) = Q 2 + Q(t)/D(t), so that for small Q, we have coD ~ 
Q{D + Q/2Q 2 ). 

Let us start by assuming our field is propagating and evolv¬ 
ing forwards (only), with perfectly matched D* fields; so that 
D = C.D 1 . but then it happens that Q changes by <5Q over 
a small interval St, likewise t; changes by St;. The D ± will 
no longer be matched, and now the total field splits into two 
parts that evolve in opposite directions. The part that contin¬ 
ues to evolve forward has D + nearly unchanged, but the for¬ 
ward evolving D has changed size (and is now — St;)) 
to stay perfectly matched according to the new Q. The rest of 
the old D (°c St;) now propagates backwards, taking with it 
a tiny fraction of the original D + (and is ; S 'q ). 

Comparing the two backward evolving D components at t 
and t + St, and taking the limit St —>• 0 enables us to estimate 
that the backward evolving D field changes according to 

2Q 

backward = ^ + ^2 [^/®] D t,forward' (B.0.2) 

Now, using the small-Q approximation for ft), we can write 

d : D QMckward = —^—2 ftQ](B.0.3) 

(ft) + ii) 




±iQD ± ± — Q. 

20 ^ 


(A.0.9) where the exponential part removes any oscillations due to the 
linear part of Q; i.e. if Q = xD then 


Appendix B: The no accumulation approximation 

In the main text, I describe the no accumulation approxima¬ 
tion in spectral terms as a RWA approximation. However, it 
is hard to set a clear, accurate criterion for the RWA approx¬ 
imation to be satisfied in the general case, since it requires 
knowledge of the entire propagation before it can be justified. 
In this appendix, I take a different approach to determine the 
conditions under which the approximation will be satisfied. 


dt D 0.backward ~ , (B.0.4) 

(ft) + Q) 

So here we see that backward evolving fields are only gener¬ 
ated from forward evolving fields due to changes in the under¬ 
lying conditions (i.e. either material response or pulse proper¬ 
ties), but that for the reflection to be strong those changes will 
have to be significant on the order of a period, or be periodic 
so that phase matching of the the backward wave could occur. 
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